Optimal evaluation of single-molecule force spectroscopy experiments 
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The forced rupture of single chemical bonds under external load is addressed. A general framework 
is put forward to optimally utilize the experimentally observed rupture force data for estimating 
the parameters of a theoretical model. As an application we explore to what extent a distinction 
between several recently proposed models is feasible on the basis of realistic experimental data sets. 
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Introduction: Single-molecule force spectroscopy [l| 
refers to the experimental observation of chemical disso- 
ciation by pulling apart the molecular complex of interest 
at a constant velocity v until the bond breaks. The eval- 
uation of the resulting rupture force data is a non trivial 
task [2, H, 0, S H, 0, Q : for one and the same w, the rup- 
ture forces are found to be randomly distributed over a 
wide range, and for different v, different such probability 
distributions are obtained. While both on the experi- 
mental and modeling sides a great amount of work has 
led to substantial progress and sophistication, much less 
effort has been spent to improve the still rather basic 
methods of connecting and comparing theory and exper- 
iment. This is the subject of our present work. 

Our starting point is the probability density v) 
that a dissociation event occurs at a pulling force /, given 
the pulling velocity v and any theoretical model with cer- 
tain model parameters /x. Now, our main question is: 
What is the optimal estimate of those model parameters 
/X that can be extracted from any given set of N rupture 
forces f = {fi}iLi and pulling velocities v = {vil'^^i ? 
Since the fi are statistically independent, the probability 
of observing the given set of rupture forces f reads 

N 

p{f\iji,w) = Ylp^{U\^l,v,) . (1) 

1=1 

The main result of our paper is that the optimal param- 
eter estimate is obtained by simply maximizing with 
respect to fi: no other "recipe" is able to yield estimates 
closer to the true parameter values systematically, i.e., 
on the average over many data sets f . 

Example: The explicit form of pi in ([l]) depends on the 
specific model one is considering, and similarly for the 
meaning and even the number of the model parameters 
fi. While our general theory applies to any model, an il- 
lustrative and particularly simple example is provided by 
the most widely used model viewing the dissoci- 

ation as a rate process of the form n{t) = —k {f{t)) n{t), 
where n{t) denotes the bond survival probability and 
k{f{t)) the dissociation rate at the instantaneouspuUing 
force f{t), and adopting the approximations P, 

fit) = Kvt , fc(/) - exp{A + af} . (2) 

Here, k is the net elasticity of the setup and kv is the so- 
called loading rate. While they are considered as known 



[H 0i M — (-^jQ;) ^re the two unknown model pa- 
rameters in the case of our specific example at hand. 
Their physical meaning is discussed in detail, e.g., in 
[l|, 0]: k{0) = exp{A} represents the force-free dissoci- 
ation rate and aksT the dissociation length (distance 
between potential well and barrier along the reaction 
pathway), where fc^T is the thermal energy. Having thus 
completely specified the model [H, 0, Q , a straightforward 
calculation yields for this particular example the explicit 
result 

Pi{I\lJ',v) = exp<^ \ . (3) 

KV y KV a ) 

Formal analysis: The quantity in ([T]) is called the like- 
lihood and plays a central role in Bayes' theorem [lo| 

p(/x|f , v) ^ p{{\fi, v) p(/x, v)/p(f , v) . (4) 

The left-hand side represents the "likeliness" of /x, given 
the data f, v, and hence is clearly of central interest 
for our purposes. Considering also the right-hand side 
as a function of /x, it is equal to the likelihood from 
^ times p(/x, v), encapsulating all our knowledge about 
/X before the measurement, times a /x-independent fac- 
tor l/p{i,v). We emphasize that we will not use the 
Bayesian formalism in our actual calculations below, only 
in their intuitive interpretation. 

Next we exploit the fact that typically a quite large set 
of rupture data f is available. Thus, focusing on large N, 
it is convenient to rewrite ([T|) as 

p(f|At,v) = exp{-7Vsw(f,/x,v)} (5) 
SAr(f,/x,v) := -N^^^\npi{fi\fM,Vi) . (6) 

i=l 

Furthermore, we assume that the relative frequency with 
which the different pulling velocities v are sampled con- 
verges toward a well-defined limit p{v) for N oo. Fi- 
nally, we assume that the rupture forces fi have been 
sampled according to the "true" distribution pi (/^ |/Xg, w^) 
with unknown, "true" model parameters /Xg. Then it fol- 
lows from the law of large numbers [ll| that 

SAr(f,/x,v) s(/x) := -(lnpi(/|/x,w)) (7) 

for iV — > oo, where {• • • ) indicates an average over / and 
V with weight pi(/|/Xq, v) p{v). Hence, s^r is an intensive, 
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entropy like quantity. Observing that s(/i) — s(/i.Q) is a 
relative entropy of the form (ln[pi(/|/XQ, i;)]), 
one can infer [111 ] that s^fi) has a unique absolute min- 
imum at = //g. For any given f and v, wc denote by 
fi* — /x*(f, v) the maximum of the likelihood p{i\fj,,v) 
with respect to fi, or, equivalently, the minimum of 
SAr(f, /x,v) in (O. Since sn converges for large to- 
ward s according to ([7]), also the minimum fi* of the 
former converges to the minimum fi^ of the latter. Con- 
sequently, for fi close to fi* and large N, we can expand 
SAr(f , fi, v) up to second order about its minimum at /x* 
and the Hessian matrix of SAr(f, /x*,v) can be replaced 
by the Hessian H — ff(/Xo) of 5(^*0)' i-^-i 

sw(f,/x* + A,v) - SAr(f,/x*,v) + . (8) 

For large this is a very good approximation for all 
/i- values with an appreciable weight in ([S]), i.e., 

p(f l/x, v) cx exp{-7V(/x - /x*)tif(/x - /x*)/2} . (9) 

Within this narrow peak region, the factor v) in 
(HI , though usually unknown in detail, can be considered 
as approximately constant, i.e., p(/x|f,v) oc p(f|/x,v). 
Given f and v, the likelihood ^ thus quantifies the "like- 
liness" that the "true" model parameters are fi. 

Upon repeating the entire set of N pulling experiments 
with the same set of pulling velocities v, a different set 
of rupture data f will be sampled, yielding a different 
maximum likelihood estimate fi*. While the probabil- 
ity distribution of f is given by ((T|) with fi = fiQ, what 
can we say about the distribution of the maximum like- 
lihood estimates fi*l To determine its first moments, we 
differentiate ^ and choose A. — fi^ — fi* , resulting in 

fj,* - fj,Q = -H-^dsNi{,fJ,o,^)/dn . (10) 

Averaging over f yields zero on the right-hand side, as 
can be inferred from ([6]), ([7]) and the fact that /1q is the 
minimum of s. Hence, 

(/^*) = /^o > (11) 

i.e., the maximum likelihood estimate is "unbiased". An 
analogous but somewhat more involved calculation [l^l 

yields for the second moments the result 

([/x*-/Xo][/x*-/Xo]t) = (iVi7)-i . (12) 

Observing that {N H)^^ is the covariance matrix of the 
distribution from ([5]), we arrive at our 

First main conclusion: For any given, sufficiently large 
data set f, the expected deviation of the concomitant 
maximum likelihood estimate fi* from the "true" param- 
eters /Xg immediately follows from the "peak width" of 
likelihood considered as a function of fi. 

Similarly, from the higher moments one can infer (l^ 
that fi* is Gaussian distributed, yielding with ^ our 

Second main conclusion: Apart from the peak posi- 
tion and a normalization factor, the likelihood ([1]) for 
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FIG. 1: Solid histogram: Numerically determined distribu- 
tion of the first components of the maxima fi* = (A*, a*) of 
the likelihood ([T]) for 10000 "computer experiments" . For each 
of them, A'^ — 400 rupture forces / were sampled according to 
Q, 100 for each of the four loading rates kv = 50, 200, 1000, 
5000 pN/s and with "true" parameters Ao = —5 and ao — 0.1 
pN~\ These are typical numbers in "real experiments" [3]. 
Thin lines: Likelihood Q for the first 15 of the 10000 exper- 
iments after integrating over a, shifting the maximum to Ao, 
and normahzing (some are almost indistinguishable). Dotted 
histogram: Distribution of the estimates for A according to 
the "standard method", as described in the main text. 

one given data set f looks practically the same as the dis- 
tribution of the maximum likelihood estimates fi* from 
many repetitions of the TV pulling experiments. 

Figure 1 illustrates these findings by means of the ex- 
ample from ([3]). Since two-dimensional distribu- 
tions are difficult to compare graphically, we focus on 
the marginal distributions for the first component A of 
fi — (A, a) (the findings for a are similar). The close 
agreement of the 15 thin lines with the histogram in Fig. 
1 very convincingly illustrates our two conclusions above. 

In view of the argument below ([5]), it seems intuitively 
quite plausible that the maximum of the likelihood fi* 
should be the best possible guess for the unknown true 
parameters /Xg. A more rigorous line of reasoning starts 
with an arbitrary "recipe" of estimating the true param- 
eters fj,Q from a given data set f , formally represented 
by some function /i(f). The only assumption is that this 
recipe is unbiased, i.e., upon repeating the same experi- 
ment many times, on the average, the "true" parameters 
are recovered, (/i(f)) = Mo- By generalizing the well- 
kown Cramer- Rao inequality llj, which in turn is basi- 
cally a descendant of the Cauchy-Schwarz inequality, one 
can show [l^ for any such "recipe" fi{f) that 

([/i-Mo][/i-Mo]t)-(Afi/)-^>0, (13) 

i.e., the matrix on the left-hand side is non- negative def- 
inite. Comparison with (jl2p yields our 

Third main conclusion: There is no unbiased estima- 
tor fi of the true parameters /.tg which on the average 
outperforms the maximum likelihood estimate fi*. 

The remaining possibility that a biased estimator may 
be even better is rather subtle to treat rigorously, but 
intuitively this seems quite unlikely. Furthermore, in the 
above conclusion we exploited the relation (I12p which is 
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FIG. 2: Rupture force distribution for different loading rates nv. Histograms: numerically generated rupture forces according 
to (|15l) with 7 = 'i-l'i, Ao = —5, ao = 0.1 pN~^, eo = 15. For each kv, we sampled 500 forces, i.e., A'^ = 2000. Solid: maximum 
likelihood fit according to (|15p for 7 = 1/2 and 7 = 2/3 (not distinguishable in this plot). Dashed: same for 7 = 1. 

Upon repeating the entire "numerical experiment", the resulting plots always look practically the same. 



strictly correct only for asymptotically large N . Finally, 
the criterion of minimizing the left-hand side in (fT3|) itself 
is in principle also debatable, but hardly in practice. Be- 
ing unable to make any further progress along these lines, 
we directly compared the maximum likelihood estimate 
with other known "recipes" of evaluating single-molecule 
rupture data 0, 0, i, H, SH, [ll . In all cases we found 
that the maximum likelihood was superior. 

Three case studies: 1. In single molecule force spec- 
troscopy, the most widely used "recipe" for estimating 
parameters consist of the following steps: (i) Fit a Gaus- 
sian to the observed rupture force distribution for a fixed 
pulling velocity v and approximate the most probable 
rupture force /* by the maximum of that Gaussian, (ii) 
Plot /* for different v versus ln(w) and fit the resulting 
points by a straight line, (iii) Assume that the model 
0, ([31) is applicable and deduce its model parameters 
fi = (A, a) from the slope and the axis intercept of the 
straight hne as detailed, e.g., in [U, 0, H H, III. We 
have applied this procedure to each of the 10000 exper- 
iments in Fig. 1 and plotted the distribution of the re- 
sulting estimates for A in Fig. 1. The conclusion is that 
the maximum likelihood estimate represents a substantial 
improvement compared to the so far "standard method" 
of data evaluation in this field. 

2. Generalizations of the rate ^ of the form 

Hf) = (1 - 7a//e)^/''"^ gA+e[i-(i-7"/A)'/-'] (14) 

with three model parameters /x = (A, a, e) have recently 
attracted considerable interest [I4j]. Here, A and a have 
the same physical meaning as in Eq. e := Ei,{0)/ ksT 
stands for the force-free activation energy barrier in units 
of the thermal energy ksT, while 7 € {1/2,2/3,1} la- 
bels three different models: For 7 = 1 the parameter e 
drops out and one recovers 7 = 2/3 reproduces the 
Kramers rate for a cubic reaction potential, and 7 = 1/2 
corresponds to a parabolic potential well with a cusp bar- 
rier |14l | . The resulting rupture force distribution 

= ^-p y-- j 

(15) 



with k{f) from (|14p can be determined analogously to 
([3|). There is an ongoing debate in the literature about 
which of the three models is most appropriate to evaluate 
experimental rupture data [l3|. Taking for granted that 
one of the three models approximates the "truth" satis- 
factorily, choosing fj. — fjb* is - according to our above 
conclusions - the closest one can get to the "full truth" 
on the basis of one given data set f . In case of disagree- 
ment about the "true" 7-value, a fully objective selec- 
tion criterion seems impossible to define in principle. In 
practice, the usual criterion is the comparison with the 
basic "true" quantity observed experimentally, namely, 
the distribution of rupture forces. In view of Fig. 2, 
we conclude that under typical experimental conditions 
it is absolutely impossible to decide whether 7 = 1/2 or 
7 = 2/3 is "better" , and even 7=1 performs almost as 
well [ll. 

3. In Fig. 3 the same comparison as in Fig. 2 is 
repeated, but now for real experimental data from d. 
Again, the models HH), ^ with 7 = 1/2 and 7 = 2/3 
are hardly distinguishable; 7 = 1 differs slightly more, 
while the "standard method" yields a completely differ- 
ent "best fit". However, none of them satisfactorily de- 
scribes the "experimental reality" . The same incompati- 
bility is recovered for all other experimental data sets we 
analyzed so far, see also Q- An almost perfect agreement 
(within the statistical uncertainty of the experimental 
data) is obtained by means of yet another recent exten- 
sion [3] of dl]), ([3]), considering the parameter a itself as 
randomly sampled from 

p{a) = (27ra2)-i/2 exp{-(a - af /2<jI} , (16) 

resulting in a model with three parameters /x = 
(A, q;,(Tc). Possible reasons for such a heterogeneity of 
the dissociation rate are uncontrollable variations of the 
experimental conditions or of the complicated biomolec- 
ular complex itself [7]. 

Conclusions: The maximum n* of the likelihood ([T]) is 
the best possible estimate for the unknown model param- 
eters /I, given an appropriate model and a (sufficiently 
large) set of rupture forces f . The accuracy of this esti- 
mate follows from the dispersion of the (approximately 
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FIG. 3: Same as Fig. 2 but for experimental rupture forces from [8|. The number of rupture events for each loading rate 
is indicated in brackets. Upper row: Experimental data (histograms) and maximum likelihood fit pi{f\fi*,v) according to 
model (I14p . (|15p for 7 = 1 (dashed), 7 = 2/3 (solid), and 7 = 1/2 (dash-dotted). Lower row: Experimental data (histograms), 
maximum likelihood fit pi(/|/i*,i;) according to model (|16|) (solid), and best fit according to the "standard method", as 
described in the main text (dotted). 



Gaussian) likelihood peak about fi*. The procedure is 
extremely simple and general. For example, the pulling 
velocities Vi may be all the same, all different, or dis- 
tributed in any other way, and the pulling force f{t) may 
or may not increase linearly with time (only in the first 
case is there a well defined loading rate k; cf. 

By means of a "least-squares fit" one gets - by defini- 
tion - the best possible agreement between theory and ex- 
periment with respect to any given "deviation measure" . 
A typical example is to optimize the agreement between 
experimental and theoretical rupture force distributions. 
In particular, the resulting agreement with the data in 
Figs. 2 and 3 would be (at least slightly) better than 
for any of the depicted theoretical lines. However, our 
present goal is not to optimally fit rupture force distri- 
butions but rather to fit the unknown model parameters 



as closely as possible: they are the quantities of prime in- 
terest, and any other kind of fitting procedure is mainly 
an intermediate step in order to estimate them. Using 
the rupture force distributions to fit parameters seems 
natural, but our paper shows that one can do better. 
Our present comparison of the rupture forces in Figs. 2 
and 3 serves a different purpose: once the parameters are 
estimated as well as possible according to our method, 
the resulting rupture force distributions can be used as 
an independent consistency test for a given hypothetical 
model. 
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